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ABSTRACT 

O 

' The Lya line transfer in an extremely thick medium of neutral hydrogen is 

' investigated by adopting an accelerating scheme in our Monte Carlo code to skip 

a large number of core or resonant scatterings. This scheme reduces computing 
^ ' time significantly with no sacrifice in the accuracy of the results. We applied 

■ this numerical method to the Lya transfer in a static, uniform, dustless, and 

plane-parallel medium. Two types of photon sources have been considered, the 
midplane source and the uniformly distributed sources. The emergent profiles 
show double peaks and absorption trough at the line-center. We compared our 
results with the analytic solutions derived by previous researchers, and confirmed 
that both solutions are in good agreement with each other. We investigated the 
directionality of the emergent Lya photons and found that limb brightening is 
observed in slightly thick media while limb darkening appears in extremely thick 
media. The behavior of the directionality is noted to follow that of the Thomson 
scattered radiation in electron clouds, because both Lya wing scattering and 
Thomson scattering share the same Rayleigh scattering phase function. The 
mean number of wing scatterings just before escape is in exact agreement with 
the prediction of the diffusion approximation. The Lya photons constituting 
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the inner part of the emergent profiles follow the relationship derived from the 
diffusion approximation. We present a brief discussion on the application of our 
results to the formation of Lya broad absorption troughs and P-Cygni type Lya 
profiles seen in the UV spectra of starburst galaxies. 

Subject headings: line: formation — radiative transfer — galaxies: starburst — 
galaxy: formation 

1. Introduction 

Recent observational achievements in spectroscopy of galaxies in the very early universe 
(Steidel et al. 1996, 1999) have enabled us to study starburst galaxies whose star formation 
rates are very high. Since these galaxies are at 3 < ^ < 5, we can observe their ultraviolet 
continua and Lya in the optical band, which can be covered by ground-based large telescopes. 
Furthermore, Lya is one of the most conspicuous emission lines, and the investigation of the 
Lya line transfer becomes important in understanding the physical nature of these galaxies. 

The Lya profiles of many starburst galaxies in the early universe are categorized into 
three types : symmetric Lya emission, asymmetric or P-Cygni type one, and broad absorp- 
tion in damping wings (Tenorio-Tagle et al. 1999). Moreover, local starburst galaxies also 
show similar characteristics (Kunth et al. 1998). The P-Cygni type Lya is beheved to orig- 
inate from the expanding neutral medium surrounding the Lya source, which is supported 
by the existence of the low ionization interstellar absorption lines blueshifted relative to the 
line-center (Kunth et al. 1998; Lee & Ahn 1998; Ahn et al. 2000). The Voigt profile fittings 
for the broad absorption lines extending to the Lorentzian wings and P-Cygni type absorp- 
tion lines revealed that their neutral hydrogen column density A^hi = 10^*^ — 10^^ cm^^ or the 
line-center optical depth tq = 10^ — 10^. However, the Voigt profile fitting is not satisfactory 
in the red part of Lya emission which is thought to be contributed by back-scattered pho- 
tons. This is exemplified by the ill-fitted red part of Lya emission in the spectrum of Haro 
2 (Legrand et al. 1997). No adequate computational method to describe the back-scattered 
photons has been developed so far. 

In extremely thick media of hydrogen, where Lya photons wander in both real space 
and frequency space, the line transfer can be described by the diffusion approximation. In 
the diffusion approximation we assume that only wing scatterings play a significant role in 
transferring Lya photons in space. However, this approximation does not fully take into 
account the fact that the photons move back and forth between core and wing wavelength 
regimes during their transfer in an extremely thick medium. While being a core photon the 
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scatterings are resonant and therefore local. Wing scatterings are associated with signifi- 
cant spatial excursion. Furthermore, the scattering phase functions for resonant scatterings 

associated with the 5*1 Pi 2. transitions are different from that associated with wing 

2 2 ' 2 

scatterings which is idential with that of the classical Rayleigh scattering. This fact should 
be considered in a very careful manner in computing the physical quantities such as the 
polarization and the angular distribution of the emergent radiation field. 

In this paper, we investigate the Lya line transfer in an extremely thick and static 
medium by using a Monte Carlo method. We will focus only on the dustless case in this 
work, even though the dusty case is more realistic. The effect of dust was briefly dealt by 
Ahn et al. (2000) and more careful investigation is deferred. Our assumption can be valid 
for the primeval objects that had not experienced any star-formation episode. 

This is the second paper in a series on the Lja profile formation in starbursting galaxies. 
In the first paper (Ahn et al. 2001) we studied the Lya transfer in moderately thick media. 
We use the Monte Carlo method, which is frequently adopted in the Lya fine transfer, because 
of its simplicity and flexibility. However, the Monte Carlo method becomes inefficient when 
the line-center optical depth tq of scattering medium is very large. We need to follow the 
scattering of each photon until the escape, and the number of scattering is proportional to tq 
for very thick case. This means that the number of computations grows in proportional to Tq 
which could vary by several orders of magnitudes depending on the situation. Our previous 
version of the code consumes unacceptably long computing time for extremely thick cases. 
We therefore extend our efforts to very large optical depth regime in this paper by developing 
the 'accelerating scheme'. 

This paper is composed as follows. In section 2, we describe our scheme to accelerate 
the Monte Carlo code. In section 3 we show our results. In section 4 we discuss possible 
applications of our code. The flnal section summarizes our major flndings. 

2. Configuration and Monte Ccirlo Method 

Since we are interested in the emergent Lya proflles in various situations, we flrst in- 
troduce the dimensionless parameter x deflned by 

X = Au/Aud = (i^ - 1^0)/ Aud, (1) 

which describes the frequency shift from the line-center vq in units of the Doppler shift 
Ai/D = J^o{vth/c)- Here Vth is the thermal speed of the scattering medium, c is the speed of 
light, and is the line-center frequency. 

The medium considered in this study is a plane-parallel slab, which is uniformly flUed 
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with pure hydrogen atoms. The optical depth at the Lya hne-center from the center to both 
sides of the medium along the normal direction is Tq. We consider two types of the source 
distribution: a source in the midplane, and sources uniformly distributed along the vertical 
line. Wc show the configuration considered in this study in Figure 1. The underlying ideas 
for these two cases is that the midplane source can model a giant H II region embedded in 
a neutral region while the uniformly distributed source can mimic the starbursting region 
which has young stars mixed with the partially ionized gas. These configurations are the 
same ones that have been considered in the previous works (Neufeld 1990; Adams 1972). 

There have been several Monte Carlo approaches to the resonance hne transfer in an 
optically thick and static medium in the literature (Avery & House 1968; Adams 1972; Meier 
& Lee 1981; Gould & Weinberg 1996). The main advantages of Monte Carlo methods lie 
in the geometrical flexibility of scattering media and the simplicity in solving the radiative 
transfer problem, compared with the direct numerical integration of the basic radiative 
transfer equation. However, one severe limitation of the technique is that a large amount 
of computing time is required for the cases with extremely high optical depths. In order 
to handle the Lya transfer for extremely large Tq using currently available computers, we 
modified the prototypical Monte Carlo code developed in our previous papers (Ahn et al. 
2000, 2001). 



2.1. Acceleration Scheme for Extremely Thick Medium 

As was shown numerically in the previous papers (Ahn et al. 2000, 2001) and explained 
in accordance with Adams (1972), in extremely thick media, Lya photons escape from the 
media by both frequency diffusion and spatial diffusion. When a Lya photon is emitted with 
the line center frequency in an optically thick medium, then it experiences a large number of 
core scatterings. During the core scatterings, the photon does not move much because the 
line-center opacity is very large. During a series of a large number of core scatterings, the 
Lya photon can be scattered by a rapidly moving atom in the tail of the Maxwcll-Boltzmann 
velocity distribution. Once the photon is scattered by such an atom, its frequency changes 
abruptly and it becomes a wing photon. For a wing photon, the hydrogen medium becomes 
more transparent and the photon traverses a longer distance. In a moderately thick medium 
with 10^ < To < 10^/a, where a is the Voigt parameter, Lya photons escape from the medium 
during at most a few wing scatterings; but in extremely thick media with otq > 10^, the 
wing photon can not escape and have to experience a number of continuous wing scatterings, 
whose number is approximately given by {N) — x\ assuming random walk approximation 
(Adams 1972). Here Xs — (otq)^/^ is the frequency at the peaks of profiles. 
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During such wing scatterings some photons become core photons again by the so called 
'restoring force' (Adams 1972), and then these photons suffer a large number of successive 
core scatterings. One cycle of these processes is called 'an excursion' (Adams 1972), and 
the photon spatially transfers by the repetition of excursions before it eventually achieves 
its escape from the medium during the wing scatterings: this final step of escape is called 
'single longest excursion'. 

Considering these transfer processes, we can see that the major consumption of time 
in the Monte Carlo code occurs in the step of the computation of the outgoing wavevector 
associated with resonant or core scatterings, whose number of successive core scattering is 
order of Tq . Therefore, the computing time can be significantly saved if we simplify this step 
in a clever way. Our code was developed to distinguish the resonant core scattering from 
the non-resonant wing scattering. Thus, in order to save the computing time, whenever the 
photon experiences a resonant scattering during its transfer, we can skip following successive 
resonant scatterings and to describe the next scattering by a fast moving hydrogen atom. 

We now explain the skipping procedure in detail. The main idea for the procedure 
is to assume that any photons with \xi\ < Xc suffer a large number of core scatterings 
locally and finally emerges from the local scattering region with frequency shift \x2\ > Xc 
after a scattering with a rapidly moving atom, where Xc is the critical frequency shift which 
defines the borderline between core and wing photons. In the accelerating scheme, we just 
skip the core scatterings and assigns X2 in an appropriate way. For the Lorentzian profile 
of Lya scattering, Xc can be safely set to < Xc < ^/7^ for moderately thick media with 
10^ < To < lO'^/a, and to Xc = \/7r for extremely thick media with aro > 10^. Here we obtain 
y/n from the Eq. (10) by substituting the effective wing optical depth = 1. We find that 
the slightly larger Xc can be permitted for even larger tq. 

In order to describe the scattering by a fast moving atom, we select the velocity (v) 
of the scattering atom from the volcano-type thermal velocity distribution (Avery & House 



where Xc < v < Vmax, rather than from a simple Gaussian distribution. Here Vmax is safely 
set to be as large as Vmax = 10. We call this method 'the accelerating scheme'. We note that 
this method is valid for otq <^ 10^ only when Xc —>■ 0, which is tantamount to restoring the 
original code. 



1968) 




(2) 
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2.2. Angular Redistribution and Polarization 

In this subsection, we provide the details of the angular redistribution and polarization 
of the scattered radiation. Our code treats the scattering of Lja photons in a manner 
very faithful to the atomic physics associated with the fine structure of hydrogen. The 
level splitting of the excited state of Lja transition, is lOGHz, which amounts to 

the Doppler width of 1.3 km s'^. Because this is much smaller than the thermal speed in 
a medium with T > 100 K, one normally neglect the ^-Pi^i level splitting, treating it as a 
single level. However, even though T < 100 K, the level splitting can be neghgible when the 
medium has a large optical thickness. 

As is described in Ahn et al. (2000, 2001), in our Monte Carlo code we compute the 
probability that the given photon is resonantly scattered with one of the two levels Pi s and 

the probability that the scattering occurs in the damping wings at each scattering event. We 
adopt the density matrix formalism to describe the angular distribution and polarization of 
the scattered Lja, where the density operator is represented by a 2 x 2 Hermitian matrix. 
The density matrix element associated with the scattered photon ppp/ is related with that of 
the incident photon paa' by 

^ Ej-EA-hw Er-EA-h. ' 

/,/'=P3/2,Pi/2,a,a' 

Here, the matrix elements due to other nP(n > 2) terms are neglected, and the radial part 
of the matrix elements < 2p\\ r ||ls > is treated as a constant. The relevant Clebsch-Gordan 
coefficients are found in Condon & Shortley (1951); Lee et al. (1994). 

As is shown by Stenflo (1980), scattering of Lja in the damping wings is characterized 
by the classical Rayleigh scattering phase function due to the quantum interference between 
the two levels Pi a. In this case, one can show that the Lya scattering in the wings can be 
regarded as a resonance scattering between two levels with J — and J = 1 as far as the 
scattering phase function is concerned. The computation of the matrix elements in this case 
is straightforward and the result is 

p'li = ^0 Pu ~ cos 9 sin 2 A0 + cos^ 9 sin^ A0 p22 

p'i2 = - cos 9' sin 2A(f) pu + (cos 9 cos 9' cos 2A0 + sin 9 sin 9' cos A0) pi2 

— cos ^^(sin sin sin A0 + - cos cos sin2A0) p22 (4) 
P22 — cos^ 9' sin^ A0 pn + cos 9' (2 sin 9 sin 9' sin A0 + cos 9 cos 9' sin 2 A0) pi2 
+(cos 9 cos 9' cos A0 + sin 9 sin 9')^^ P22, 
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where the incident radiation is characterized by the wavevector kj = (sin 6 cos 0, sin 9 sin 0, cos 9) 
and the outgoing wavevector k/ is correspondingly given with angles 9' and (j)' with A0 = 
0' — 0. Here, we assume that the circular polarization is always zero due to the azimuthal 
symmetry in the scattering geometry and initial condition. 

In the case of resonance scattering between J — \ and J — \ levels, a straightforward 
calculation shows that the density matrix elements associated with the outgoing photons are 
given by 

p'li = (5 + 3cos2A0)pii 

+[(5-3 cos 2A0) cos^ 6* + 2 sin^ 9]p22 - 6 cos 6' sin 2A0pi2 
p'i2 — 3 sin 2A0 cos 9' pn + 6(cos 9 cos 9' cos 2A0 + sin 9 sin 9' cos A0)pi2 

+3 cos 9{—2 sin 9 sin 9' sin A0 — cos 9 cos 9' sin 2A0)p22 

= [(5-3cos2A0)cos2^' + 2sin^^')Pii 

+ [(5 + 3 cos 2A0) cos^ 9 cos^ ^' + 2 cos^ 9 sin^ ^' 

+ 12 cos A0 cos 6*' cos 9 sin 6* sin 6*' + 2 cos^ 9' sin^ 6* + 8 sin^ 9 sin^ 6'']p22 

+ (6 sin 2 A0 cos 9 cos^ 6*' + 2 sin A0 cos 6*' sin 6* sin 9') p^- (5) 

When the scattering is resonant with the level Pi, then the scattered radiation is charac- 
terized with the isotropic angular distribution and hence completely unpolarized. Therefore, 
the density matrix is constant and | times the identiy matrix irrespective of the incident 
and outgoing wavevectors. 

The density matrix elements associated with the outgoing photon given in Eqs. (3-5) 
are normalized by the unit trace condition. It is readily verified that a 90° scattering of 
a completely unpolarized incident photon yields 100 percent polarization in the case of a 
wing scattering and the degree of polarization in the case of a resonant scattering between 
Si Ps becomes | (e.g. Chandrasekhar (1960); Lee et al. (1994)). 

In Figure 2 we compare the results obtained by using the accelerated code with those 
calculated by using the original code in our previous works. Here we set the temperature 
of the scattering media to be T = 10 K or a = 1.47 x 10~^, and the photon source is 
located at the origin of the central plane of the slab, which is the same configuration to the 
midplane sources by symmetry. In the figure the solid lines stand for those results without 
the accelerating scheme, and the dashed lines represent the results using the modified code 
with Xc = 1. Notable deviation between the two results is seen for arg < 1.47 x 10^ or 
O'To < 2.6 X 10^. Here tq is the line-center optical depth defined in terms of the normalized 
Voigt function in the previous researches (Adams 1972; Neufeld 1990), where Tq = x/ttto. 
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The original code can be regarded as a special case of the accelerated code when 0. 
In order to check this, we calculate the emergent profile for the smaller core-wing boundary 
frequency, Xc- For the case with Tq = 10^, we show our result with Xc = 0.6 by the dotted 
line in Figure 2. Note that the dashed line for Tq = 10"^ is obtained for x^ = 1. In the figure 
we confirm that the emergent profile converges to that from the original code as Xc 0. It 
is clear that the accelerated code is much more efficient than the original code, and it gives 
the correct results. 



3. Line Transfer in Extremely Thick Media 

3.1. Emergent Profiles 

In this subsection, we compare our results with analytic solutions and check the validity 
of our accelerating scheme. In extremely thick media with otq > 10^, wing scatterings 
dominate the line transfer in space and the Voigt function H{x, a) can be approximated by 

Hix,a)^-^^, (6) 

where the optical depth at a frequency shift x is given by — TQH[a,x). Adams (1972); 
Harrington (1973); Neufeld (1990) introduced the diffusion approximation according to which 
the Lya photons transfer in space only by wing scatterings, and the distance traversed 
per scattering is described by Eq. (6). Neufeld (1990) derived an analytic solution for the 
extremely thick media, 

\/6 x^ 1 

J(±To, X) = — — ^^gj^ j^^4 1^^^ 1/2 ( I ^3 1 /q^^)] • (^) 

We perform Monte Carlo calculations for a monochromatic source with its input fre- 
quency = 0. We simulate several cases with different aro values so that both moderately 
thick and extremely thick cases are included in the computations. In each case, the source 
is located at the midplane of a static neutral slab. 

First, we set the temperature of the scattering medium to be T = 10 K, which corre- 
sponds to a = 1.49 X 10^^. The relevant situation for this case may be found in the later 
evolutionary epoch of starburst galaxies, when the cold supcrshcU with large H 1 column 
density expands very slowly. The column density (A^hi) of normal galaxies and dwarf galax- 
ies ranges Ahi = 10^^ — 10^^ cm~^, which corresponds to the Lya hne-center optical depth 
To = 10^ — 10^. The line-center optical depth is related to the H I column density A^hi via 

— 1/11 rTi-l/2 r A^l 

To = 1.41 



mi 



10i3cm-2 



(8) 
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where T4 = T/IO^K and T is the temperature of the scattering medium. 

Here we choose Tq = 10^, 10^ 10^ 10^, 10^ or aro = 1.49 x 10, 1.49 x 10^, 1.49 x 10^ 
1.49 X 10^, 1.49 X 10^. We note that these are transformed respectively into afo = 2.64 x 
10, 2.64 X 10^, 2.64 x 10^ and 2.64 x lOl 

Our results are shown in Figure 3a and Figure 3b. We see that the profiles of our 
Monte Carlo results are in good agreement with the analytic solutions for aro > 10^. For 
To = 10^, we also show the emergent profile for the case of = 0.6 which is closer to the 
result of original code in Ahn et al. (2000, 2001). We see that it also shows deviation from 
the analytic solution derived by Neufcld (1990), which is expected due to the invalidity of 
diffusion approximation for ar^ < 10^. We conclude that our accelerating scheme is valid for 
the various range of optical depths. 

Next, we perform Monte Carlo calculations for hotter media, whose Voigt parameter is 
given by a = 4.71 x 10~^ corresponding to T = K. Wc choose tq = 10^, 10'', 10^, and 
10^, which correspond to aro = 4.71 x 10^4.71 x 10^,4.71 x 10*^, and 4.71 x 10^ We note 
that these are transformed respectively into afo = 8.35 x 10^, 8.35 x 10^, 8.35 x 10"^, and 
8.35 X 10^. In Figure 4, we show the results of our accelerated Monte Carlo code and compare 
them with the analytic solution, Eq. (7). The solid lines represent the results of our Monte 
Carlo calculations, and the dotted fines do analytic solutions. We see that our numerical 
results are in very good agreement with analytic solutions for very thick media (tq > 10^). 
In extremely thick media, we can safely assume that only wing scatterings contribute to the 
spatial transfer of Lya photons, and so our accelerated code can describe the line transfer 
process very well. 

Third, we also perform calculations for sources uniformly distributed along the vertical 
direction of the slab. The relevant appfication for this case can be found in the giant H II 

region, where the Lya sources and the scattering media may be mixed together. We assume 
that the typical temperature of H II region is T = 10*^ K, and that the size of H II region 
is kpc scale. In this environment the degree of ionization can lead to the volume density of 
neutral hydrogen nni = 1 cm~^. Hence the total column density of the H II region can be 
A'hi ~ 10^° cm~^, which corresponds to the line-center optical depth of ro ~ 10"^ — 10^. In the 
code, we generate the sources at the location {px,Py,Pz) = (0,0, r^), where r^ = {2R — l)ro 
and R is the uniform random number having range of [0, 1]. 

When photon sources arc located near the surface, the optical depth toward the nearest 
boundary is much less than that toward the midplane of the slab. Hence, we adjust the core- 
wing boundary frequency {xc) in accordance with the vertical location of photon sources 
which is denoted by Pz'- we used Xc — A{1 — Pz/tq), where A is a core- wing boundary 
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frequency in units of thermal velocity for midplane sources. This formula gives correct 
behaviors for both thin and thick cases. Clearly — > as jp^j — > Tq. 

In Figure 5 we show our results only for the extremely thick cases with tq = 10^, 10^ and 
a — 4.71 X 10~^. We note that the results for the moderately thick cases can be found in the 
previous paper (Ahn et al. 2001). Our results (sohd lines) are compared with the analytic 
solutions (dashed lines). From the analytic solutions for the case of midplane source (Neufeld 
1990), we obtain the analytic solutions for the case of distributed sources by convolving the 
midplane solution, Eq. (7), with the uniform source positions (ts), 

/TO 
Jm{±To,X;Ts)S{Ts)dTs. (9) 
-To 

In the figure we can see that our numerical results agree very well with the analytic ones, 
and that the peaks are more shifted to the line-center compared to the midplane cases. This 
is because the sources near the surface contribute to the central part of emergent profiles. 

Adams (1972) argued that in extremely thick media with aro > 10^, the emergent flux 
depends on the product aro, not on the individual parameters. This is evidently seen in 
Eq. (7), in which J{±Tq,x) is dependent upon arg. In order to check this fact, we perform 
Monte Carlo calculations for a few pairs of a and tq with common aro values. For aro = 
4.71 X 10^ (or afo = 8.35 x 10^), we perform calculations for two pairs; (a = 4.71 x 10~^, 
To = 10^) and (a = 4.71 x 10"^, To = 10^). For otq = 4.71 x 10=^ (or afo = 8.35 x 10^), 
we choose (a = 4.71 x 10-^ To = 10^) and (a = 4.71 x 10"^ To = 10^). The other pair is 
(a = 4.71 X 10-^ To = 10^) and (a = 4.71 x 10"^, tq = 10^), which have aro = 4.71 x 10^. 
Our results are shown in Figure 6, where we see that, for the extremely thick cases aro > 10^, 
the line profiles for each pair with the same aro agree very well with each other. Thus we 
confirm that our code reproduces the results in Adams (1972) very well. 

Finally, we also consider the radiative transfer of contimium photons near Lya. We 
simply assume that an input spectrum is fiat. In Figure 7 we show our results for tq = 
10^, 10^, 10^, 10'' fixing a = 1.49 x 10~^, where symmetric double peaks and one absorption 
trough at the line-center are seen. The absorption trough broadens as Tq gets larger. We 
show the wavelength (A) in units of A on the top edge of the figure, which is scaled by 
A oc T^/^. Here T is the temperature of the scattering medium. Note again that the spectra 
are symmetric relative to the line-center. 
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3.2. Limb brightening and darkening 



In this subsection we consider the directionahty of emergent Lja photons. In the Monte 
Carlo code we collect emerging photons with the directional information into the equally 
spaced bins in = cos^^ where 9 is the angle between the normal direction of the slab 
and the wavevector of the photons. Since the number of photons along is proportional 
to I{ijL)cos9dQ, where Q is the sohd angle, we obtain the intensity /(//) by dividing the 
number of photons by We call the function /(//) the 'directionahty'. This definition of 
directionality is different from that used by Phillips & Meszaros (1986, hereafter PM86), 
who defined it as the intensity normalized with respect to the intensity normal to the slab. 

The results are shown in Figure 8. The curves for aro < 10^ are nearly isotropic, because 
most Lya photons escape from the medium by single longest flights: each Lya photon is 

scattered locally near the source and becames a wing photon by frequency diffusion, when 
the medium gets transparent. Therefore, a large number of local scatterings isotropize the 
radiation field, which is characterized by no or (at least a negligible) polarization. 

As To > 10^, the directionality gradually converges to that of the Thomson-scattered ra- 
diation emergent from a Thomson-thick electron medium represented by the thick solid line. 
This Thomson limit was obtained by Chandrasekhar (1960) and was confirmed numerically 
by PM86. PM86 also provided the results for optically thin as well as optically thick cases. 
Therefore, we compare our results with those from the radiative transfer in a pure electron 
scattering medium, because Thomson scattering shares the same Rayleigh phase function as 
Lya wing scattering. 

It is convenient to consider the effective wing optical depth, r^, which is the optical 
depth at the effective frequency of emergent photons, Xs = (aro)^^^. From Eq. (6) we can 
see that 



We compare this optical depth with the optical depth for the Thomson scattering defined 
by Te = CTNg, where ctt = 0.665 x 10"^"^ cm^ is the Thomson scattering cross section, and 
Ne is the electron column density. 

When aTQ is small, after a large number of local resonant scatterings, a Lya photon 
becomes a wing photons due to the scattering with a fast moving hydrogen atom. Now, for 
this wing photon the slab is a transparent medium, especially along the normal direction 
and escape to this direction occurs without subsequent scatterings. However, to the grazing 
direction, because of simple geometrical considerations, there can be a significant probability 
of further wing scatterings before escape. Therefore, the situation becomes quite similar to 
the electron scattering in a plane-parallel slab that is quite thin and illuminated from the 




(10) 
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other side. Limb brightening is a natural consequence because we see more sources in the 
grazing direction, as is shown in Figure 8. 

For the extremely thick cases with aro > 10^, the effective wing optical depth is esti- 
mated to be > 6 by Eq. (10). This corresponds to the result of the Thomson scattered 
radiation in an electron cloud with its optical depth Te > 6. This limiting behavior of the 
directionality may be summarized as 



which is an approximate relation proposed by PM86 for the numerical solution introduced in 
the Table XXIV of Chandrasekhar (1960). This is shown by a thick solid line in Figure 8. We 
see that the directionality of Lya, being transferred in an extremely thick hydrogen medium, 
converges to that for the case of optically thick electron cloud. In this optical depth regime, 
the optical depth for Lya wing photons becomes large, and the photons spatially transfer 
in the medium by random walks. Hence, the escape of photons occurs preferentially to the 
normal to the slab, where the opacity is relatively small. We regard this phenomenon as a 
limb darkening, which was also called a 'beaming' by PM86. 

In an extremely thick medium, Lya photons experience a number of successive wing 
scatterings just before their escape from the medium during their single longest excursion. 
How many times do the photons experience the successive wing scatterings just before es- 
cape? According to Adams (1972), in an extremely thick medium with an isotropic Lya 
source located at the midplane, photons acquire a typical frequency shift Xs = {aroY^^ dur- 
ing excursions (see also section 2 of this paper). In an extremely thick medium, the spatial 
transfer of Lya photons occurs during a series of wing scatterings, and so we can approxi- 
mate the spatial transfer by random walks. Thus from Eq. (10) we can estimate the mean 
number of successive wing scatterings just before escape {{Nyj)) by 



In the accelerated Monte Carlo calculations we count the number of wing scatterings 
just before escape, and compare the result with the above relationship. In the upper panel 
of Figure 9 we show the spectral distribution of the number of wing scatterings just be- 
fore escape. We can sec that the photons emerging with larger frequency shift have been 
experienced more wing scatterings just before escape than those with smaller frequency shift. 




(11) 




(12) 



or 




(13) 



-13- 



In order to check the relationship of Eq. (13) , we show in Figure 10 the graph of log^g (N^) 
vs. log^o (^^0, where (N^) is the value averaged over frequency, and aro represents the impor- 
tance of wing scatterings in Lya line transfer. Here the number of photons for each point is 
4000, and the error in A^^ is the flux weighted Poisson error. The least square fitting gives 
us the following relationship, 

logio(A^«,) = (0.674 ± 0.010) logio otq - (0.552 ± 0.042). (14) 

We sec that the power and the coefficient in Eqs. (13) and Eq. (14) are in good agreement 
with each other. 



4. Discussion 

We have investigated the Lya line transfer in extremely thick, dustless, uniform, and 
static media. In order to reduce the calculation speed, we developed an accelerating scheme 
for our previous Monte Carlo code, in which wc skip an extremely large number of core 
scatterings that may occur during the line transfer. The code also deals with the angular 
redistribution accurately by using the density matrix formalism where the phase functions 
calculated in a manner faithful to the atomic physics. The modified code passes the require- 
ment in both efficiency and accuracy when calculating the Lya line profiles emergent from 
extremely thick media. We checked that the line profile is dependent only upon aro, which 
was argued by Adams (1972). With the code we obtained solutions for both the midplane 
source and the uniformly distributed sources, both of which show excellent agreement with 
the analytic solutions derived by Neufeld (1990). In the latter case, the peaks tend to shift 
toward the line-center due to the sources located near the slab surface. We also presented our 
results on the radiative transfer of continuum photons near Lya, and the emergent spectra 
show the double peaks and the broad absorption at the line-center. However, it is natural to 
think that dust extinction will destroy a significant fraction of photons forming the double 
peaks, and eventually a broad absorption trough will form. 

We also examined the directionality of emergent Lya photons. When the medium is 
slightly thick with otq << 10^, we observe a hmb brightening effect, because the slab becomes 
transparent to wing photons. As aro increase, the directionality eventually converges to that 
of the Thomson scattered radiation in a thick electron cloud which was investigated by 
Chandrasekhar (1960). The limiting behavior of the curves is the beaming of Lya photons 
into the direction normal to the slab. In short, the Lya emission from the star-forming 
regions that are surrounded by scattering neutral media shows the limb brightening when 
utq -C 10"^, and the limb darkening when otq > 10^. These facts can be physically explained. 
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since both the hya scattering by neutral hydrogen at wings and the Thomson scattering in 
an electron cloud have the same Rayleigh phase function. 

We also counted the mean number of successive wing scatterings just before escape, 
which is denoted by (Nyj), and by Monte Carlo calculations we confirmed the relationship 
(Nw) — ^{O'ToY^^, which is derived from the diffusion approximation. 

In conclusion, it is important to consider the quantum mechanical properties of Lya 
scatterings and the exact transfer processes, and our accelerated Monte Carlo code can 
afford those requirements as well as the computing time. 

Our works can be applied to local and distant starburst galaxies. Approximately half 
of their Lya lines show broad absorption troughs in the Lorentzian wings. These absorption 
troughs are believed to be formed by optically thick, static H I media with some amount of 
dust. However, the media in this study have been assumed to be dustless, which is unrealistic 
for almost all situations. Hence, we will take dust effects into account in the future works. 
Moreover, the emergent profiles are also influenced by the physical quantities such as porosity 
of dust and homogeneity of the medium as well as dust abundance and H I column density, 
so that wc are now undertaking investigations about the effects of these components on the 
line formation. 

Another half of the starburst galaxies, either nearby or remote, show P-Cygni type Lya 
emission. As Kunth et al. (1998) suggested, expanding media around the star-forming regions 
play an important role in avoiding the destruction of Lya photons by dust and forming the 
P-Cygni type Lya profiles. When we calculated the Lya line formation in an expanding 
medium, it is difficult to deal with the back-scattered photons. We have been investigating 
this topic by applying the formalism developed in this paper, and the results will be reported 
in another paper. 

In conclusion, Lya from starburst galaxies bears important and useful information on 
the physical environment of those systems including the kinematics, H I column density, 
abundance and porosity of dust in the surrounding media of the star-forming regions in those 
galaxies. Being very luminous in addition to being sensitive to many physical quantities, 
Lya can be easily observed by using the modern observational technology. 

This work was achieved as a part of doctorial dissertation in the School of Earth and 
Environmental Sciences of Seoul National University financially supported by Brain Korea 
21 of the Korean Ministry of Education. SHA and HML appreciate the support. Also this 
work is completed with the financial support of the Korea Institute for Advanced Study, 
which is pleasantly acknowledged. 
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Fig. 1. — Configuration studied in this work. We consider the static and uniform slab. We 
also consider two types of source distribution: in one type the Lya source is located at the 
origin denoted in the figure by O, and in the other they are uniformly scattered along z-axis. 
The vertical thickness of the slab is 2to, and we define // as the cosine of the angle between 
the wave vector of the emergent photon and the slab normal. 
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Fig. 2. — Comparisons between the results of the original code (solid line) and those of 
the code with the accelerating scheme (dashed line). The bottom horizontal axis represents 
frequency shift x in units of the thermal width, and the top horizontal axis represents wave- 
length for T — lOK. Comparisons were made for the case of midplane sources radiating 
monochromatic Lya photons with xq — 0. The dotted line represents the result for a small 
core- wing boundary frequency, Xc- 
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Fig. 3. — Lya profiles emerging from thick media of different tq's, where we fix a = 1.49 x 10"^ 
or T = 10 K. The bottom horizontal axis represents frequency in units of the thermal 

width, and the top horizontal axis represents wavelength shift for T = lOK. Here the 
solid lines stand for results of our Monte Carlo calculations, and the dashed lines for the 
analytic solutions given by Neufeld (1990). The total flux of the line is normalized to I/Att 
in accordance with Neufeld's normalization. We note that the profiles are symmetric about 
the origin, x = 0. 
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Fig. 4. — Lycc profiles emerging from extremely thick media with different Tq's and a common 
Voigt parameter a — 1.49 x 10~^ corresponding to T = 10 K. Here the solid lines stand for 
the results of our Monte Carlo calculations, and the dashed fines for the analytic solutions 
given by Neufeld (1990). Here horizontal axis represents frequency shift x in units of the 
thermal Dopplcr width, and the total flux of the line is normalized to l/47r in accordance 
with Neufeld's normalization. Note that the profiles are symmetric about the origin, x — Q. 
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Fig. 5. — Our emergent profiles (solid lines) are compared with Neufeld's analytic solution 
(dotted lines). 
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Fig. 6. — Emergent profiles for the cases of distributed sources. Comparing with the mid- 
plane solutions of the same a and Tq (dotted lines), their peaks are shifted to the line-center. 
Here solid fines represent our Monte Carlo results, and the dashed lines do the analytic 
solutions obtained by convolving Neufeld's solution. 
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Fig. 7. — Emergent profiles of a few pairs of Monte Carlo simulations which have the common 
utq values. The pair of physical parameters in the calculations are listed in the text. We 
see that the cases with the same otq show the same emergent profiles. Note that the profile 
is symmetric with respect to the origin x = 0. On the top horizontal axis is shown the 
wavelength (A) in units of A, which scales A oc T^/^. Here T is the temperature of the 
scattering medium, and set to be T = 10"^K or a = 4.71 x 10^''. For T = lOK or a = 1.49 x 
10^^, the wavelengh scales are smaller by a factor of vTOOO = 31.6, because Avd oc T^/^. 




Fig. 8. — Emergent spectra near Lyct when the continuum source is at the center of the 
scattering medium. Here a = 1.49 x 10~^ is assumed, and the spectra are normahzed to the 
continuum leveL We also denote the wavelength in A on the top edge of the figure. 
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Fig. 9. — Directionality of emergent Lyo; photons for various optical depths Tq and Voigt 
parameters a of scattering media. We define the directionality by the flux divided by The 
solid thick line represents the limiting behavior of directionality for the Thomson scattered 
radiation in an opaque electron cloud with the radiation source deep in the cloud. The 
radiation distribution is almost isotropic when the wing scattering optical depth is small, 
and converges to take the form of a linear function J(/x) = a + fc/x, which is also the limit 
attained in the very thick Thomson scattering medium obtained by Chandrasekhar (1960). 
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Fig. 10. — Numbers of successive wing scatterings just before escape of Lya photons, Nyj 
versus the frequency, x. The source is located in the midplane of the slab whose temperature 
is T = lOK or a = 1.49 x 10"^. The line-center optical depths are tq = 10^ 10^ 10^. We 
show A^^ with the Poisson errors, and for clarity we added one and two to the curves for 
To = 10^ and tq = 10^, respectively. 




Fig. 11. — The mimbcr of a scries of wing scatterings just before escape {Ny^) versus effective 
optical tliickness (aro). Both are presented in logarithmic scales. The error bars are those 
of flux- weighted errors, and the straight line shows the least square fit. 



